﻿#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#pragma warning(disable:4996)
#include <stdbool.h>
#include<math.h>
#include <float.h>

#define MAX_POINTS 1000  // 假设最多1000个数据点
#define MOMENTUM 0.9
#define MAX_ITERATIONS 10000000
#define LEARNING_RATE 0.00001  // 更小的学习率
#define CONVERGENCE_THRESHOLD 1e-8

typedef struct {
    float sse;      // 剩余平方和
    float r_squared; // 决定系数
    double chi_square;
    double adj_r_squared;
} FitGoodness;

typedef struct {
    float a;       // 斜率
    float b;        // 截距
    float c;
    float d;
    float sse;     // 最终剩余平方和
    int iterations;// 实际迭代次数
} LinearFitResult;

void show(float list[], int length)
{
    for (int i = 0; i < length; i++)
    {
        printf("%.2f\n ", list[i]);
    }
    printf("\n");
}

int read_data_from_file(const char* filename, float x[], float y[]) {
    FILE* file = fopen(filename, "r");
    if (file == NULL) {
        printf("无法打开文件: %s\n", filename);
        return -1;
    }

    int count = 0;
    char line[256];  // 足够存储一行数据

    while (fgets(line, sizeof(line), file) != NULL && count < MAX_POINTS) {
        // 更灵活的格式匹配，处理可能的空格
        if (sscanf(line, " (%f , %f )", &x[count], &y[count]) == 2 ||
            sscanf(line, "(%f,%f)", &x[count], &y[count]) == 2 ||
            sscanf(line, "%f %f", &x[count], &y[count]) == 2) {
            count++;
        }
        else {
            // 尝试更宽松的解析方式
            char* ptr = line;
            char* end;

            // 解析第一个浮点数
            x[count] = strtof(ptr, &end);
            if (ptr == end) continue; // 转换失败

            // 跳过分隔符
            ptr = end;
            while (*ptr && (*ptr == ',' || *ptr == ' ' || *ptr == '\t')) ptr++;

            // 解析第二个浮点数
            y[count] = strtof(ptr, &end);
            if (ptr == end) continue; // 转换失败

            count++;
        }
    }

    fclose(file);

    if (count == 0) {
        printf("警告: 文件 %s 中没有找到有效数据\n", filename);
    }
    else {
        //printf("成功从 %s 读取 %d 组数据\n", filename, count);
    }

    return count;
}

int leastSquareLinearFit(float x[], float y[],
    int num, float* a, float* b, FitGoodness* goodness)
{
    if (num < 2) {
        printf("至少需要2个数据点\n");
        return -1;
    }

    float sum_x = 0, sum_y = 0, sum_xy = 0, sum_x2 = 0;
    float y_mean = 0;

    // 计算各项累加和及y的均值
    for (int i = 0; i < num; i++)
    {
        sum_x += x[i];
        sum_y += y[i];
        sum_xy += x[i] * y[i];
        sum_x2 += x[i] * x[i];
        y_mean += y[i];
    }
    y_mean /= num;

    float denom = num * sum_x2 - sum_x * sum_x;
    if (denom == 0)
    {
        printf("无法计算，分母为零\n");
        return -1;
    }

    // 计算回归系数
    *a = (num * sum_xy - sum_x * sum_y) / denom;
    *b = (sum_x2 * sum_y - sum_x * sum_xy) / denom;

    // 计算拟合优度指标
    if (goodness != NULL)
    {
        float sse = 0.0f;   // 剩余平方和
        float tss = 0.0f;   // 总平方和

        for (int i = 0; i < num; i++)
        {
            float y_pred = (*a) * x[i] + (*b);  // 预测值
            sse += powf(y[i] - y_pred, 2);      // 残差平方
            tss += powf(y[i] - y_mean, 2);      // 总平方
        }

        goodness->sse = sse;
        goodness->r_squared = 1.0f - (sse / tss);
    }

    return 0;
}

bool LeastSquares(float* x, float* y, unsigned int data_num,
    float* a, float* b, float* c, FitGoodness* goodness)
{
    // 检查数据点数量
    if (data_num < 3) {
        printf("至少需要3个数据点进行二次拟合\n");
        return false;
    }

    // 计算各项累加和
    double sumx = 0, sumx2 = 0, sumx3 = 0, sumx4 = 0;
    double sumy = 0, sumxy = 0, sumx2y = 0;
    double y_mean = 0;

    for (int i = 0; i < data_num; i++) {
        double xi = x[i];
        double xi2 = xi * xi;
        double xi3 = xi2 * xi;
        double xi4 = xi2 * xi2;

        sumx += xi;
        sumx2 += xi2;
        sumx3 += xi3;
        sumx4 += xi4;
        sumy += y[i];
        sumxy += xi * y[i];
        sumx2y += xi2 * y[i];
        y_mean += y[i];
    }
    y_mean /= data_num;

    // 构建正规方程矩阵
    double matrix[3][4] = {
        {data_num, sumx, sumx2, sumy},
        {sumx, sumx2, sumx3, sumxy},
        {sumx2, sumx3, sumx4, sumx2y}
    };

    // 高斯消元法解方程组
    for (int k = 0; k < 3; k++) {
        // 部分主元选择
        int max_row = k;
        for (int i = k + 1; i < 3; i++) {
            if (fabs(matrix[i][k]) > fabs(matrix[max_row][k])) {
                max_row = i;
            }
        }

        // 交换行
        if (max_row != k) {
            for (int j = k; j < 4; j++) {
                double temp = matrix[k][j];
                matrix[k][j] = matrix[max_row][j];
                matrix[max_row][j] = temp;
            }
        }

        // 消元
        for (int i = k + 1; i < 3; i++) {
            double factor = matrix[i][k] / matrix[k][k];
            for (int j = k; j < 4; j++) {
                matrix[i][j] -= factor * matrix[k][j];
            }
        }
    }

    // 回代求解系数
    *c = matrix[2][3] / matrix[2][2];
    *b = (matrix[1][3] - matrix[1][2] * (*c)) / matrix[1][1];
    *a = (matrix[0][3] - matrix[0][1] * (*b) - matrix[0][2] * (*c)) / matrix[0][0];

    // 计算拟合优度指标
    if (goodness != NULL) {
        double sse = 0.0;    // 剩余平方和
        double tss = 0.0;    // 总平方和
        double chi_sq = 0.0; // 卡方统计量

        for (int i = 0; i < data_num; i++) {
            double xi = x[i];
            double y_pred = (*c) * xi * xi + (*b) * xi + (*a); // 预测值
            double residual = y[i] - y_pred; // 残差

            sse += residual * residual; // 残差平方和
            tss += (y[i] - y_mean) * (y[i] - y_mean); // 总平方和

            // 卡方统计量（假设误差方差为1）
            if (fabs(y_pred) > 1e-10) { // 避免除以零
                chi_sq += (residual * residual) / fabs(y_pred);
            }
        }

        goodness->sse = sse;
        goodness->r_squared = 1.0 - (sse / tss);
        goodness->chi_square = chi_sq;
    }

    return true;
}

bool LeastSquaresCubic(float* x, float* y, unsigned int n,
    float* a, float* b, float* c, float* d,
    FitGoodness* goodness)
{
    if (n < 4) {
        printf("至少需要4个数据点进行三次拟合\n");
        return false;
    }

    // 计算各项累加和
    double sum_x = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0, sum_x5 = 0, sum_x6 = 0;
    double sum_y = 0, sum_xy = 0, sum_x2y = 0, sum_x3y = 0;
    double y_mean = 0;

    for (int i = 0; i < n; i++) {
        double xi = x[i];
        double xi2 = xi * xi;
        double xi3 = xi2 * xi;
        double xi4 = xi2 * xi2;
        double xi5 = xi3 * xi2;
        double xi6 = xi3 * xi3;

        sum_x += xi;
        sum_x2 += xi2;
        sum_x3 += xi3;
        sum_x4 += xi4;
        sum_x5 += xi5;
        sum_x6 += xi6;
        sum_y += y[i];
        sum_xy += xi * y[i];
        sum_x2y += xi2 * y[i];
        sum_x3y += xi3 * y[i];
        y_mean += y[i];
    }
    y_mean /= n;

    // 构建正规方程矩阵
    double A[4][5] = {
        {n,      sum_x,  sum_x2, sum_x3, sum_y},
        {sum_x,  sum_x2, sum_x3, sum_x4, sum_xy},
        {sum_x2, sum_x3, sum_x4, sum_x5, sum_x2y},
        {sum_x3, sum_x4, sum_x5, sum_x6, sum_x3y}
    };

    // 高斯消元法解方程组
    for (int k = 0; k < 4; k++) {
        // 部分主元选择
        int max_row = k;
        for (int i = k + 1; i < 4; i++) {
            if (fabs(A[i][k]) > fabs(A[max_row][k])) {
                max_row = i;
            }
        }

        // 交换行
        if (max_row != k) {
            for (int j = k; j < 5; j++) {
                double temp = A[k][j];
                A[k][j] = A[max_row][j];
                A[max_row][j] = temp;
            }
        }

        // 消元
        for (int i = k + 1; i < 4; i++) {
            double factor = A[i][k] / A[k][k];
            for (int j = k; j < 5; j++) {
                A[i][j] -= factor * A[k][j];
            }
        }
    }

    // 回代求解系数
    *d = A[3][4] / A[3][3];
    *c = (A[2][4] - A[2][3] * (*d)) / A[2][2];
    *b = (A[1][4] - A[1][2] * (*c) - A[1][3] * (*d)) / A[1][1];
    *a = (A[0][4] - A[0][1] * (*b) - A[0][2] * (*c) - A[0][3] * (*d)) / A[0][0];

    // 计算拟合优度指标
    if (goodness != NULL) {
        double sse = 0.0;    // 剩余平方和
        double tss = 0.0;    // 总平方和
        double chi_sq = 0.0; // 卡方统计量

        for (int i = 0; i < n; i++) {
            double xi = x[i];
            double y_pred = (*d) * xi * xi * xi + (*c) * xi * xi + (*b) * xi + (*a); // 预测值
            double residual = y[i] - y_pred; // 残差

            sse += residual * residual; // 残差平方和
            tss += (y[i] - y_mean) * (y[i] - y_mean); // 总平方和

            // 卡方统计量（假设误差方差为1）
            if (fabs(y_pred) > 1e-10) { // 避免除以零
                chi_sq += (residual * residual) / fabs(y_pred);
            }
        }

        goodness->sse = sse;
        goodness->r_squared = 1.0 - (sse / tss);
        goodness->chi_square = chi_sq;
        goodness->adj_r_squared = 1.0 - (1.0 - goodness->r_squared) * (n - 1) / (n - 4);
    }

    return true;
}

void gradientDescentLinearFit(float* x, float* y, int n, LinearFitResult* result) {
    if (n < 2) {
        printf("至少需要2个数据点\n");
        return;
    }

    // 特征缩放：计算x的均值和范围
    float x_mean = 0.0f, x_min = x[0], x_max = x[0];
    for (int i = 0; i < n; i++) {
        x_mean += x[i];
        if (x[i] < x_min) x_min = x[i];
        if (x[i] > x_max) x_max = x[i];
    }
    x_mean /= n;
    float x_range = x_max - x_min;

    // 初始化参数
    float a = 0.0f, b = 0.0f;
    int iter = 0;

    for (iter = 0; iter < MAX_ITERATIONS; iter++) {
        float grad_a = 0.0f, grad_b = 0.0f;

        // 计算梯度（使用缩放后的x值）
        for (int i = 0; i < n; i++) {
            float x_scaled = (x[i] - x_mean) / x_range; // 特征缩放
            float error = (a * x_scaled + b) - y[i];
            grad_a += error * x_scaled;
            grad_b += error;
        }

        // 更新参数
        float new_a = a - LEARNING_RATE * grad_a / n;
        float new_b = b - LEARNING_RATE * grad_b / n;

        // 检查收敛
        if (fabs(new_a - a) < CONVERGENCE_THRESHOLD &&
            fabs(new_b - b) < CONVERGENCE_THRESHOLD) {
            a = new_a;
            b = new_b;
            break;
        }

        a = new_a;
        b = new_b;
    }

    // 将参数转换回原始尺度
    a = a / x_range;
    b = b - a * x_mean;

    // 存储结果
    result->a = a;
    result->b = b;
    result->iterations = iter + 1;

    // 计算最终SSE
    result->sse = 0.0f;
    for (int i = 0; i < n; i++) {
        float error = (a * x[i] + b) - y[i];
        result->sse += error * error;
    }
}

void GradientDescentQuadraticFit(float* x, float* y, int n, LinearFitResult* result) {
    if (n < 3) {
        printf("至少需要3个数据点\n");
        return;
    }

    // 特征缩放：计算均值和标准差
    float x_mean = 0.0f, x_sq_mean = 0.0f;
    for (int i = 0; i < n; i++) {
        x_mean += x[i];
        x_sq_mean += x[i] * x[i];
    }
    x_mean /= n;
    x_sq_mean /= n;
    float x_std = sqrtf(x_sq_mean - x_mean * x_mean);
    if (x_std < 1e-10) x_std = 1.0f;

    // 初始化参数和动量
    float a = 0.0f, b = 0.0f, c = 0.0f;
    float ma = 0.0f, mb = 0.0f, mc = 0.0f; // 动量项
    int iter = 0;

    for (iter = 0; iter < MAX_ITERATIONS; iter++) {
        float grad_a = 0.0f, grad_b = 0.0f, grad_c = 0.0f;
        float max_grad = 0.0f;

        // 计算梯度
        for (int i = 0; i < n; i++) {
            float x_scaled = (x[i] - x_mean) / x_std;
            float x_sq = x_scaled * x_scaled;
            float error = (a * x_sq + b * x_scaled + c) - y[i];

            grad_a += error * x_sq;
            grad_b += error * x_scaled;
            grad_c += error;

            // 跟踪最大梯度分量用于自适应学习率
            max_grad = fmaxf(max_grad, fabsf(error * x_sq));
            max_grad = fmaxf(max_grad, fabsf(error * x_scaled));
            max_grad = fmaxf(max_grad, fabsf(error));
        }

        // 自适应学习率
        float learning_rate = LEARNING_RATE;
        if (max_grad > 0) {
            learning_rate = LEARNING_RATE / max_grad;
        }

        // 应用动量更新
        ma = MOMENTUM * ma + (1 - MOMENTUM) * grad_a / n;
        mb = MOMENTUM * mb + (1 - MOMENTUM) * grad_b / n;
        mc = MOMENTUM * mc + (1 - MOMENTUM) * grad_c / n;

        // 更新参数
        float new_a = a - learning_rate * ma;
        float new_b = b - learning_rate * mb;
        float new_c = c - learning_rate * mc;

        // 检查收敛
        if (iter > 1000 &&
            fabs(new_a - a) < CONVERGENCE_THRESHOLD &&
            fabs(new_b - b) < CONVERGENCE_THRESHOLD &&
            fabs(new_c - c) < CONVERGENCE_THRESHOLD) {
            a = new_a;
            b = new_b;
            c = new_c;
            break;
        }

        a = new_a;
        b = new_b;
        c = new_c;
    }

    // 将参数转换回原始尺度
    float x_mean_sq = x_mean * x_mean;
    result->a = a / (x_std * x_std);
    result->b = (b - 2 * a * x_mean / x_std) / x_std;
    result->c = c - b * x_mean / x_std + a * x_mean_sq / (x_std * x_std);
    result->iterations = iter + 1;

    // 计算SSE
    result->sse = 0.0f;
    for (int i = 0; i < n; i++) {
        float error = (result->a * x[i] * x[i] + result->b * x[i] + result->c) - y[i];
        result->sse += error * error;
    }
}

void gradientDescentCubicFit(float* x, float* y, int n, LinearFitResult* result) {
    if (n < 4) {
        printf("至少需要4个数据点\n");
        return;
    }

    // 特征缩放：计算均值和标准差
    float x_mean = 0.0f, x2_mean = 0.0f;
    for (int i = 0; i < n; i++) {
        x_mean += x[i];
        x2_mean += x[i] * x[i];
    }
    x_mean /= n;
    x2_mean /= n;
    float x_std = sqrtf(x2_mean - x_mean * x_mean);
    if (x_std < 1e-10) x_std = 1.0f;

    // 初始化参数和动量项
    float a = 0.0f, b = 0.0f, c = 0.0f, d = 0.0f;
    float ma = 0.0f, mb = 0.0f, mc = 0.0f, md = 0.0f;
    int iter = 0;
    float best_sse = FLT_MAX;
    float best_a = a, best_b = b, best_c = c, best_d = d;

    for (iter = 0; iter < MAX_ITERATIONS; iter++) {
        float grad_a = 0.0f, grad_b = 0.0f, grad_c = 0.0f, grad_d = 0.0f;
        float max_grad = 0.0f;
        float current_sse = 0.0f;

        // 计算梯度和SSE
        for (int i = 0; i < n; i++) {
            float x_scaled = (x[i] - x_mean) / x_std;
            float x2 = x_scaled * x_scaled;
            float x3 = x2 * x_scaled;
            float error = (a * x3 + b * x2 + c * x_scaled + d) - y[i];

            grad_a += error * x3;
            grad_b += error * x2;
            grad_c += error * x_scaled;
            grad_d += error;

            current_sse += error * error;

            // 跟踪最大梯度分量
            max_grad = fmaxf(max_grad, fabsf(error * x3));
            max_grad = fmaxf(max_grad, fabsf(error * x2));
            max_grad = fmaxf(max_grad, fabsf(error * x_scaled));
            max_grad = fmaxf(max_grad, fabsf(error));
        }

        // 保存最佳参数
        if (current_sse < best_sse) {
            best_sse = current_sse;
            best_a = a;
            best_b = b;
            best_c = c;
            best_d = d;
        }

        // 自适应学习率
        float learning_rate = LEARNING_RATE;
        if (max_grad > 0) {
            learning_rate = LEARNING_RATE / max_grad;
        }

        // 应用动量更新
        ma = MOMENTUM * ma + (1 - MOMENTUM) * grad_a / n;
        mb = MOMENTUM * mb + (1 - MOMENTUM) * grad_b / n;
        mc = MOMENTUM * mc + (1 - MOMENTUM) * grad_c / n;
        md = MOMENTUM * md + (1 - MOMENTUM) * grad_d / n;

        // 更新参数
        float new_a = a - learning_rate * ma;
        float new_b = b - learning_rate * mb;
        float new_c = c - learning_rate * mc;
        float new_d = d - learning_rate * md;

        // 检查收敛
        if (iter > 1000 &&
            fabs(new_a - a) < CONVERGENCE_THRESHOLD &&
            fabs(new_b - b) < CONVERGENCE_THRESHOLD &&
            fabs(new_c - c) < CONVERGENCE_THRESHOLD &&
            fabs(new_d - d) < CONVERGENCE_THRESHOLD) {
            a = new_a;
            b = new_b;
            c = new_c;
            d = new_d;
            break;
        }

        a = new_a;
        b = new_b;
        c = new_c;
        d = new_d;
    }

    // 使用最佳参数（防止过冲）
    a = best_a;
    b = best_b;
    c = best_c;
    d = best_d;

    // 将参数转换回原始尺度
    float x_mean2 = x_mean * x_mean;
    float x_mean3 = x_mean2 * x_mean;
    result->a = a / (x_std * x_std * x_std);
    result->b = (b - 3 * a * x_mean / x_std) / (x_std * x_std);
    result->c = (c - 2 * b * x_mean / x_std + 3 * a * x_mean2 / (x_std * x_std)) / x_std;
    result->d = d - c * x_mean / x_std + b * x_mean2 / (x_std * x_std) - a * x_mean3 / (x_std * x_std * x_std);

    result->iterations = iter + 1;

    // 计算最终SSE
    result->sse = 0.0f;
    for (int i = 0; i < n; i++) {
        float error = (result->a * x[i] * x[i] * x[i] +
            result->b * x[i] * x[i] +
            result->c * x[i] +
            result->d) - y[i];
        result->sse += error * error;
    }
}

int main(int argc, char* argv[]) {
    if (argc < 2) {
        printf("用法:\n");
        printf("  gendat 输出文件名 多项式次数 系数1 系数2 ...\n");
        printf("  fitting 数据文件名 方法(最小二乘法1s/梯度下降法gd) 多项式次数\n");
        return 1;
    }

    if (strcmp(argv[1], "gendat") == 0) {
        if (argc < 4) {
            printf("gendat 用法: gendat 输出文件名 多项式次数 系数1 系数2 ...\n");
            return 1;
        }

        FILE* file1 = fopen(argv[2], "w");
        if (file1 == NULL) {
            printf("无法打开文件 %s\n", argv[2]);
            return 1;
        }

        int times = atoi(argv[3]);
        if (times < 1 || times > 3) {
            printf("多项式次数必须是1, 2或3\n");
            fclose(file1);
            return 1;
        }

        if (argc < 4 + times) {
            printf("参数不足，需要 %d 个系数\n", times + 1);
            fclose(file1);
            return 1;
        }

        float data[4] = { 0 };
        for (int i = 0; i <= times; i++) {
            data[i] = atof(argv[4 + i]);
        }

        if (times == 1) {
            for (float x = -100; x < 101; x++) {
                fprintf(file1, "(%.2f,%.2f)\n", x, data[0] * x + data[1]);
            }
        }
        else if (times == 2) {
            for (float x = -100; x < 101; x++) {
                fprintf(file1, "(%.2f,%.2f)\n", x, data[0] * x * x + data[1] * x + data[2]);
            }
        }
        else if (times == 3) {
            for (float x = -100; x < 101; x++) {
                fprintf(file1, "(%.2f,%.2f)\n", x, data[0] * x * x * x + data[1] * x * x + data[2] * x + data[3]);
            }
        }
        printf("数据已保存到 %s\n", argv[2]);
        fclose(file1);
    }
    else if (strcmp(argv[1], "fitting") == 0) {
        if (argc < 5) {
            printf("fitting 用法: fitting 数据文件名 方法(1s/gd) 多项式次数\n");
            return 1;
        }

        float x[MAX_POINTS], y[MAX_POINTS];
        int count = read_data_from_file(argv[2], x, y);
        if (count <= 0) {
            return 1;
        }

        FitGoodness goodness;
        int t = atoi(argv[4]);
        if (t < 1 || t > 3) {
            printf("多项式次数必须是1, 2或3\n");
            return 1;
        }

        if (strcmp(argv[3], "1s") == 0) {
            float a, b, c, d;
            if (t == 1) {
                leastSquareLinearFit(x, y, count, &a, &b, &goodness);
                printf("\ny= %.2f x+ %.2f, count = %d\n", a, b, count);
                printf("拟合优度:\n");
                printf("  剩余平方和(SSE): %.2f\n", goodness.sse);
                printf("  决定系数(R²): %.2f\n", goodness.r_squared);
            }
            else if (t == 2) {
                LeastSquares(x, y, count, &a, &b, &c, &goodness);
                printf("\ny= %.2f x^2+ %.2f x+%.2f, count = %d\n", c, b, a, count);
                printf("拟合优度指标:\n");
                printf("  剩余平方和(SSE): %.4f\n", goodness.sse);
                printf("  决定系数(R²): %.4f\n", goodness.r_squared);
                printf("  卡方统计量(χ²): %.4f\n", goodness.chi_square);
            }
            else if (t == 3) {
                LeastSquaresCubic(x, y, count, &a, &b, &c, &d, &goodness);
                printf("\ny= %.2f x^3+ %.2f x^2+%.2f x+%.2f, count = %d\n", d, c, b, a, count);
                printf("拟合优度指标:\n");
                printf("  剩余平方和(SSE): %.6f\n", goodness.sse);
                printf("  决定系数(R²): %.6f\n", goodness.r_squared);
                printf("  调整决定系数: %.6f\n", goodness.adj_r_squared);
                printf("  卡方统计量(χ²): %.6f\n", goodness.chi_square);
            }
        }
        else if (strcmp(argv[3], "gd") == 0) {
            LinearFitResult result;
            if (t == 1) {
                gradientDescentLinearFit(x, y, count, &result);
                printf("拟合结果: y = %.2fx + %.2f\n", result.a, result.b);
                printf("迭代次数: %d\n", result.iterations);
                printf("剩余平方和: %.6f\n", result.sse);
                // 计算R²
                float y_mean = 0.0f;
                for (int i = 0; i < count; i++) y_mean += y[i];
                y_mean /= count;

                float tss = 0.0f;
                for (int i = 0; i < count; i++) tss += powf(y[i] - y_mean, 2);

                printf("决定系数R²: %.6f\n", 1.0f - (result.sse / tss));
            }
            else if (t == 2) {
                GradientDescentQuadraticFit(x, y, count, &result);
                printf("拟合结果: y = %.4fx^2 + %.4fx + %.4f\n", result.a, result.b, result.c);
                printf("迭代次数: %d\n", result.iterations);
                printf("剩余平方和: %.6f\n", result.sse);
                // 计算R²
                float y_mean = 0.0f;
                for (int i = 0; i < count; i++) y_mean += y[i];
                y_mean /= count;

                float tss = 0.0f;
                for (int i = 0; i < count; i++) tss += powf(y[i] - y_mean, 2);

                printf("决定系数R²: %.6f\n", 1.0f - (result.sse / tss));
            }
            else if (t == 3) {
                gradientDescentCubicFit(x, y, count, &result);
                printf("拟合结果: y = %.4fx³ + %.4fx² + %.4fx + %.4f\n",
                    result.a, result.b, result.c, result.d);
                printf("迭代次数: %d\n", result.iterations);
                printf("剩余平方和: %.6f\n", result.sse);
                float y_mean = 0.0f;
                for (int i = 0; i < count; i++) y_mean += y[i];
                y_mean /= count;

                float tss = 0.0f;
                for (int i = 0; i < count; i++) tss += powf(y[i] - y_mean, 2);

                printf("决定系数R²: %.6f\n", 1.0f - (result.sse / tss));
            }
        }
        else {
            printf("未知的方法: %s (必须是1s或gd)\n", argv[3]);
            return 1;
        }
    }
    else {
        printf("未知命令: %s (必须是gendat或fitting)\n", argv[1]);
        return 1;
    }

    return 0;
}